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Abstract 

A three-dimensional numerical solver based on finite-difference solution of three-dimensional elasto- 
dynamic equations in generalized curvilinear coordinates has been developed and used to generate 
data such as radial and tangential stresses over various gear component geometries under rotation. 
The geometries considered are an annulus, a thin annular disk and a thin solid disk. The solution is 
based on first principles and does not involve lumped parameter or distributed parameter systems 
approach. The elastodynamic equations in the velocity-stress formulation that are considered here 
have been used in the solution of problems of geophysics where non-rotating Cartesian grids are 
considered. For arbitrary rotating geometries, these equations along with the appropriate bound- 
ary conditions have been cast in generalized curvilinear coordinates 1 in the present study. The 
generalized curvilinear coordinate grids for the geometries considered here are created using a new 
and enhanced elliptic grid generation algorithm 2 . The solver, Finite-Difference Development of Lin- 
ear Elasticity (FiDDLE), has been validated by comparing the predictions with the corresponding 
closed form axisymmetric steady state solutions for these rotating geometries. The presence of the 
modelled Coriolis forces in the numerical solution has been demonstrated to be appreciable at high 
rotational speeds. The theoretical solutions are, therefore, in deficit by this Coriolis effect at high 
rotational speeds since the theory excludes its presence. These high rotational speeds are present 
in, e.g., turbo-jet engines and, therefore, the Coriolis forces produced in the turbine blades need to 
be modelled. 


Keywords: 

Rotating gear component geometries, Three-dimensional elastodynamic equations, Stress tensor trans- 
formation, Generalized Curvilinear Coordinates, Finite- difference solution 


1 Introduction 

The new elliptic grid generation methodology 1 ’ 2 makes it possible to generate two-dimensional and 
three-dimensional grids automatically around or inside arbitrary geometries, without any need for 
human intervention, as has been the case till now. For dynamically changing shapes, grids can be 
regenerated automatically during simulation as and when required. This would be pertinent to cases 
where a given geometry undergoes deformation such as bending and twist. In the present application 


1 U.S. provisional patent application has been filed 
2 U.S. Patent application has been filed 



of rotating gear component geometries, gear teeth could be subject to deformation, pitting, wear, 
and eventually cracks, and grids would need to be regenerated due to the corresponding changes 
in the boundary configuration. Having resolved the difficulty of automatically updating the grids 
during a simulation, the finite-difference modeling of such problems has now become attractive. 

The need for numerical simulation of dynamic stresses over gears in mesh in both normal and 
damaged states has been delineated by the lack of any normal vibration data or any anomalous 
vibration data that may reflect the presence of gear damage such as in a pinion of a OH- 58 helicopter 
transmission. Such simulated “clean” and “fault data” will aid in developing and enhancing fault- 
detection algorithms for such rotating systems by first, calibrating the fault-detection algorithms 
with the simulated clean vibration data and second, by validating these detection algorithms against 
simulated data corresponding to known anomalies. In these simulations, various anomalies can 
be seeded and allowed to propagate in time and the corresponding data generated. In general, 
in such simulations, a vast variety of faults can be introduced into the system of interest, and 
the corresponding data could be recorded for use in the developmental work on fault-detection 
algorithms. 

Toward this end, as a first step, the present three-dimensional solver has been developed and 
validated by comparing its predictions with the corresponding theoretical steady state solutions. 
The agreement is shown to be good. In an ongoing study to be reported shortly, the structural 
dynamic solver, FiDDLE, will be validated by comparing its predictions with known theoretical 
time-dependent solutions. 


2 Governing Equations 

The three-dimensional linear elastodynamic equations of motion describing the principle of momen- 
tum conservation and the constitutive equations governing the wave phenomena within an isotropic 
elastic body can be written as a system of nine equations, three for the velocities and six for the 
stresses 3,4,6 , respectively. The velocity equations are given by 

pdtQi = djTij + fi ( 1 ) 

where the velocity vector, q t = (u, v, w). The body force vector is given by fi = (f x , f y , f z ), and 
the symmetric stress tensor, r,j, has six distinct components. The stress tensor is expressed by the 
following tensorial equation. 


dtTij = A SijdivQ + ii(djqi + dtfj) (2) 

where div Q is the divergence of the velocity vector, Q ( q z ), 5ij is the Kronecker delta, and where 
A and p are the Lame constants - these elastic constants (e.g., p = r xy le xv is the rigidity or the 
elastic shear modulus, e xy is the strain caused by the stress r xy in the xy plane) characterizing the 
elastic behavior of the body are related to the Young’s modulus of elasticity E, the Poisson’s ratio 
<7 and the bulk modulus k by the following relations. 

E = p(2p + 3A)/ (p + A) 
ct = A/ (2 (p + A)) 
k = E/ (3(1 - 2a)) 


The elastic quantities p and A are functions of space for a nonhomogeneous body. 
In the rotating frame of reference, the velocity equations (1) become 6,7 


pdtqi — &j7~ij T fi; ( | fi x v | * /2) “H qj uj -(- fi ( 3 ) 

where |fi| is a constant rotational speed (fi = Wfc), r = (x,y,z) is a positional vector, second and 
third terms on the right hand side represent the centrifugal and Coriolis forces associated with the 
rotating frame of reference. 

In generalized orthogonal curvilinear coordinates, Equations (2) and (3), can be shown to assume 
the following flux-conservative form: 
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( 4 ) 


d T (A/J) 

+dd(A/J)tit - (1 + ZyC + £ Z D)] 

+dn[(A/J)rit - (1/J)(t] x B + r\ y C + r] z D)\ 

+dd(A/JXt - (l/J)«xB + CyC + C Z D)] = R/J 

where 

A — [pU. pV ; P^b T X x-) T~xy: Ucz; Tyyi Tyz: Pzz] 

B — T X y, T xz: (A + 2//)u, //U, pW, A U, 0, Au] 

C = [r X3/ ,T w ,T yz , Au,pw,0, (X + 2p)v,pw, Xv] T 

D = [T xz ,T yz ,T zz ,\w,0,p,u,\w,p,v,(^ + 2p)w] T 

Jacobian of the transformation, J, is given by d(x,y,z)/d(t;,r],(), and the metric quantities, £ t , 
£ x , etc. have their usual meanings. The right hand side column vector, R, contains the centrifugal 
and Coriolis terms and additionally terms containing spatial derivatives of the elastic constants, for 
nonhomogeneous bodies. 

Physical quantities are normalized with the Young’s modulus, E, the acoustic speed, \J (E/p), 
and the characteristic dimension such as the radius, r, for a solid disk or a shaft and ( r ou t — ri n ) for 
an annulus or an annular disk. 


3 Boundary Conditions 

The two-dimensional and three-dimensional validation examples presented here have corresponding 
axisymmetric steady state theoretical solutions 8,9 , and the predictions are directly compared with 
these theoretical solutions for a rotating annulus (two-dimensional), rotating thin annular disk (three- 
dimensional), rotating thin solid disk (three-dimensional) and a rotating long solid shaft (three- 
dimensional). Boundary conditions corresponding to these closed form solutions are transformed 
from the generalized coordinate space, (£,p, C), to Cartesian coordinates, (x,y,z), using contravariant 
tensor transformation. For example, second order stress tensor transformation between (x,y,z) space 
and (£, p, () space is 


T ij 


EE 


d X i dx j _ik 

dx l dx k 


where \ l = (£, p, C)> b j span (£, p,£) space and 1, k span (x,y,z) space. 

The contravariant stress tensor transformation given above does not include normalizing factors. 
Using normalization, this transformation in matrix form can be written as 


Mx = r 

where the metric coefficient matrix, 
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solution vector, 

^ [Uca:; T X y : T xz: Tyy : Ty Z: T zz \ 

and the right-hand side vector, 

r = [|V£| 2 r a , |V£||Vp|r £7) , |V£||VC|t £C , |V p| 2 r w , |Vp||V<|t„ c , |VC| 2 t c< ] T 
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In a two-dimensional Cartesian coordinate system (y,z), the stresses would be related to their 
counterparts in the generalized coordinates ( rj , () as 

T vv = J 72 [|Vf?| 2 C z 2 t vv - 2\Vr]\\V C\VzCzT vi + |V<| ^Acd 

T yz = jTjt - CyCzTrjri + I \ I VCI (Vy Cz + r tzCy) T r)i ~ I^CI r ty T !z T C,c\ 

and 

T ZZ = -^2 [| Vr)\ 2 (y 2 T vv - 2|V»?||VCK<W + |VC| 2 % 2 tcc] 

where the Jacobian of the transformation, J', is given by = VyCz ~ VzCyi and V is the 

gradient operator. Here, the orthogonal generalized coordinates, (r],C), correspond to the polar 
coordinates (0, r). 

In general, for a three-dimensional case, the stress tensors in the two coordinate systems would 
be related by the matrix form given above, i.e., 


x = M x r 

The velocity boundary conditions are either of the Dirichlet type or of the Neumann type; the 
latter being derived from the governing equations for the velocity vector, once the stress tensor is 
updated at the boundaries according to the preceding formulation. 


4 Numerical Method 

A second-order in time and space, time-staggered leap frog method is used to integrate the velocity 
and stress equations. Both the spatial derivatives and time derivatives are discretized using central 
differences. A very small numerical damping term is used to eliminate the mesh drifting (checker- 
board) instability. Since we are interested in sampling rates of the order of 50 Khz for the flight 
data of interest, the time step restriction imposed by the CFL condition of hyperbolic systems akin 
to explicit methods for integration of the dynamic system are of the same order as those posed by 
the flight data sampling rate. For example, for a steel shaft of radius equal to 20 cm. rotating at 
100 revolutions per second (rps), a typical time step would be of the order of a microsecond, which 
translates to a simulation sampling rate on the order of 100 Khz. 


5 Results 

Some steady state theoretical results from linear elasticity were used as validation and verification 
base to compare the computational results with. Validation was done using three test cases: a 
rotating annulus (two-dimensional) , a rotating thin annular disk (three-dimensional) and a rotating 
thin solid disk (three-dimensional). 


5.1 Rotating Annulus 

The thickness of the axisymmetric rotating annulus is assumed to be sufficiently small as compared 
to its radius so that the radial and tangential stresses do not vary over its thickness. This case is 
thus referred to as a rotating annulus here to avoid any confusion that may arise when we discuss the 
three-dimensional rotating annular disk with small finite thickness where the radial and tangential 
stresses do vary over the thickness of the disk. The closed form solution for the rotating annulus is 
given in Ref. 8. The weight of the annulus is neglected. 

Results for the annulus shown in Fig. 1, with inner and outer radii of 10 cm. and 30 cm. 
respectively, are shown for various rotational speeds in Figs. 2 and 3. Rotational speeds of 2000 
rps, 100 rps and 5 rps are considered to identify the effect of Coriolis forces that become appreciable 
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Figure 1: A representative 72x21 grid of an annulus: 72 grid points in the circumferential direction 
and 21 points in the radial direction. 

at high rotational speeds. Fig. 2(a) shows a comparison between the analytical and the computed 
results at 100 rps. The abscissa shows the coordinate along any radial line from the inner radius to 
the outer radius, and the ordinate shows the normalized stress. The analytical radial and tangential 
stress distributions along this radial ray are shown in solid and dotted lines respectively, and the 
computed radial and tangential stress distributions are shown in dashed and dotted lines with crosses 
respectively. Agreement between the computed and the analytical results is satisfactory. The grid 
used for this case is 72x11, 11 points in the radial direction and 72 points along the circumferential 
direction. As the grid is refined to 72x21, an improvement is seen in the predictions, and the 
comparison with the theoretical results is good, as is shown in Fig. 2(b). 

Results for the rotational speed of 5 rps with a grid of 72x21 are shown in Fig. 3(a). Agreement 
between the predictions and the theoretical results is good, just as in the case of 100 rps. Next, 
results for a high rotational speed of 2,000 rps with the same grid (72x21) are shown in Fig. 3(b). 
While agreement between the predicted and theoretical radial stress distributions continues to be 
good, as for lower rotational speeds, the tangential stress distributions show a discrepancy. This 
discrepancy is clearly due to the presence of Coriolis forces that are modelled in the computations 
and thus reflected in the predictions, but which have not been accounted for in the theory. To 
confirm the presence of this Coriolis effect, another computation was made at 2000 rps, with the 
Coriolis terms turned off, and this discrepancy disappeared and the overall agreement for tangential 
stresses was as good as at lower rotational speeds. 

5.2 Rotating Thin Annular Disk 

An approximate theoretical solution of a rotating thin annular disk is given in Ref. 9. The problem 
is treated as that of plane stress, in which the only nontrivial stress components are the radial and 
the tangential stresses as in the case of the two-dimensional annulus, but now they also have a weak 
dependence on the thickness of the disk (along the axial direction). The thickness of the disk is 
taken to be 2.5 cm., while the inner and outer radii of the disk are 10 cm. and 30 cm. respectively, 
just as in the case of the two-dimensional annulus. The theoretical solution is in defect near the 
ends, but is good at axial stations removed from the end planes. Here, the radial stress does not 
vanish along the outer or inner radii, as in the case of the annulus, but the resultant radial tension 
between any two planes not too close to the end planes along the inner and outer radii vanishes. 

The grid used for the disk, shown in Fig. 4, is 72x21x7, with 7 points along the axial direction, 
21 points along the radial and 72 points along the circumferential direction. 

Results at three axial stations are shown in Fig. 5. Fig. 5(a) shows the stresses at the plane 
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(a) 72x11 grid 


(b) 72x21 grid 


Figure 2: Annulus at 100 rps: comparison between predicted and theoretical axisymmetric stress 
distributions; radial stresses given by lower set of curves, theoretical stresses in lines and predictions 
in lines with crosses 




(a) 5 rps (b) 2000 rps 


Figure 3: Annulus, 72x21 grid: comparison between predicted and theoretical axisymmetric stress 
distributions; radial stresses given by lower set of curves, theoretical stresses in lines and predictions 
in lines with crosses 
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Figure 4: A representative 72x21x7 grid of a thin annular disk: 72 grid points in the circumferential 
direction, 21 points in the radial direction and 7 points in the axial direction. 


next to the end plane; fig. 5(b) shows the stresses at the section quarter length removed from the 
end plane; fig. 5(c) shows the stresses at the mid-section of the disk. As before, predictions are 
shown with lines with crosses. The rotational speed of the disk is taken as 100 rps, which is of the 
same order as the rpm rate of the pinion gear for a OH-58 helicopter transmission. The agreement 
between the predictions and the theory is quite good. 


5.3 Rotating Thin Solid Disk 

Again, an approximate theoretical solution of a rotating solid thin disk is also given in Refs. 8 and 
9. This case is similar to the thin annular disk with the difference that there is no inner radius and 
therefore the attendant radial tension vanishes at the inner radius. 

A non-orthogonal cross-sectional grid is generated for the solid disk to avoid a polar coordinate 
singularity at the origin. Since the governing equations used in this study are strictly valid for 
an orthogonal curvilinear coordinate system, there is bound to be some discrepancy between the 
predictions and theory. The choice of the non-orthogonal cross-sectional grid in this case will thus 
help quantify the solution errors associated with departure from a strictly orthogonal grid. It should 
be noted that a strictly orthogonal grid can also be generated for the solid disk cross-section, but 
it will involve a zonal grid approach, which does not form the subject of this study. However, 
this problem can be resolved by re-casting the governing equations in non-orthogonal curvilinear 
coordinate systems, a study to be reported later. Since we are mainly interested in annular grids for 
gear studies, the coordinate singuarity problem is absent, and hence single-zone strictly orthogonal 
grids can be generated. 

The first grid, a coarse grid used for the thin solid disk is taken to be 21x21x7, with 7 points 
along the axial direction, 21 points along the radial and circumferential directions. The radius of 
the disk is considered to be 10 cm. and its thickness 2.5 cm. The grids are progressively refined 
from 21x21x7 to 41x41x11 to 61x61x11 to 81x81x11 to see the improvement in the predictions. 

Stresses at the last axial station (end plane) are shown in Fig. 6. Fig. 6(a) shows the radial 
stresses and fig. 6(b) shows the tangential stresses. Boundary conditions at the end planes are 
imposed from theory. Stresses are shown with lines with points and symbols for numerical boundary 
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Figure 5: Thin annular disk at 100 rps, 72x21x7 grid: comparison between predicted and theoretical 
axi-symmetric stress distributions; radial stresses given by lower set of curves, theoretical stresses in 
lines and predictions in lines with points 
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Radial Coordinate 


(a) radial stress 



Radial Coordinate 


(b) tangential stress 


Figure 6: Thin solid disk at 100 rps; axi-symmetric stress distributions at the end plane - boundary 
conditions imposed from theory at the end planes 


condition and the theoretical stresses are shown in solid line. The rotational speed of the disk is 
taken as 100 rps. 

Fig. 7 shows the predicted stresses at the quarter-length plane. Predictions are shown with 
lines with points and symbols and the theoretical stresses in solid line. Agreement between the 
predictions and the theoretical results improves as the grid is refined from 21x21x7 to 41x41x11 to 
61x61x11 to 81x81x11. At 81x81x11 grid resolution, the predictions agree very well with the theory. 
As anticipated, any errors in the predictions at this stage may be due to the grid nonorthogonality, 
as discussed earlier, but it should be noted that the theoretical solution itself is also approximate. 

Fig. 8 shows the corresponding results at the mid-plane for the various grids. As expected, the 
predictions are progressively improved as the grid is refined. 

Fig. 9a shows the cross-section of a 61x61x11 grid for the thin solid disk, and Fig. 9b shows a 
view of this three-dimensional grid. 
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(a) radial stress (b) tangential stress 

Figure 7: Thin solid disk at 100 rps; comparison between predicted and theoretical axi-symmetric 
stress distributions at the quarter-length plane - predictions are shown for grid resolutions, 21x21x7, 
41x41x11, 61x61x11 and 81x81x11; theoretical stress shown in solid line 



(a) radial stress 



(b) tangential stress 


Figure 8: Thin solid disk at 100 rps; comparison between predicted and theoretical axi-symmetric 
stress distributions at the quarter-length plane - predictions are shown for grid resolutions, 21x21x7, 
41x41x11, 61x61x11 and 81x81x11; theoretical stress shown in solid line 
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Figure 9: A representative 61x61x11 grid of a thin solid disk: 61x61 grid for the cross-section at any 
axial station and 11 points along the axial direction. 

5.4 Concluding Remarks 

A three-dimensional finite-difference elastodynamics simulator in generalized curvilinear coordinates 
has been developed and validated by comparing the predictions with steady-state axisymmetrical 
analytical solutions. The simulator, FiDDLE (Finite Difference Development of Linear Elasticity in 
Generalized Coordinates), will be used to generate time-dependent stress data corresponding to a 
variety of boundary conditions on a rotating gear simulating a gear pair in mesh. These numerically 
generated vibrations data will, in turn, be used to calibrate and validate fault-detection algorithms 
to be developed for fault detection in rotating systems such as helicopter and aircraft transmissions. 
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formulation that are considered here have been used in the solution of problems of geophysics where non-rotat- 
ing Cartesian grids are considered. For arbitrary geometries, these equations along with the appropriate boundary 
conditions have been cast in generalized curvilinear coordinates in the present study. 
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